##########################################################################################################################################
###### This file estimates and saves the betas based on daily data
###############################

# set to your working directory:
setwd("C:\\Users\\Admin\\OneDrive\\Desktop\\Market Power and Systematic Risk - Data and Code")
getwd()

###########################
# these are files from CRSP that contain the date (yyyymmdd), PERMNO, and holding period return
# - for each stock
Daily_Data = read.table("output.txt", header=T, sep="\t", dec=".")
# - and for the Market
Mkt_Daily = read.table("Mkt_Daily.txt", header=T, sep="\t", dec=".")
######
######
# these files are downloaded from Kenneth French's data library
Daily_FF = read.table("F-F_Research_Data_Factors_daily.txt", sep="", header=T, dec=".")


library(rmatio)

#Retain only the risk-free rate:
Daily_FF = Daily_FF[,c(1,5)]
Daily_FF$RF = Daily_FF$RF/100
#Calculate Betas (HIST) and other Betas inside of this function
#Indicate the estimation window by the amount of months:

#save Sec_ID:
Sec_ID = colnames(Daily_Data)[-1]
#remove Sec_ID
Daily_Data = Daily_Data[-1,]

Tage = 365

#Make sure Daily_Data, Daily_FF and Mkt_Daily has the same number of rows!!

RF = Daily_FF$RF[Daily_FF$Date%in%Daily_Data[,1]]
RF = RF[Daily_FF$Date[Daily_FF$Date%in%Daily_Data[,1]]%in%Mkt_Daily$DATE]

Mkt_Daily = Mkt_Daily[Mkt_Daily$DATE%in%Daily_Data[,1],]
Daily_Data = Daily_Data[Daily_Data[,1]%in%Mkt_Daily$DATE,]
#Do excess returns:
Mkt_Daily$vwretd  = Mkt_Daily$vwretd - RF
Mkt_Daily$vwretdx = Mkt_Daily$vwretx - RF

Excess <- function(x){
  Temp <- x - RF
  return(Temp)
}

#calculate excess returns (daily)
Daily_Data[,2:ncol(Daily_Data)] <- apply(Daily_Data[,2:ncol(Daily_Data)], 2, Excess)


df1 <- data.frame(Date = as.Date(as.character(Daily_Data[,1]),format="%Y%m%d"))


####################################################### Optimized Solution #####################################################################


#Make sure Daily Data and Mkt_Daily are both of the same length

Daily_Data = Daily_Data[as.Date(as.character(Daily_Data[,1]),format="%Y%m%d") %in%  as.Date(as.character(Mkt_Daily$DATE),format="%Y%m%d"),]
Mkt_Daily = Mkt_Daily[as.Date(as.character(Mkt_Daily$DATE),format="%Y%m%d") %in% as.Date(as.character(Daily_Data[,1]),format="%Y%m%d"),]

dim(Daily_Data)
dim(Mkt_Daily)
library(data.table)
Daily_Data$vwret = Mkt_Daily$vwretd

setDT(Daily_Data)

library(roll)
cols <- colnames(Daily_Data)[2:(ncol(Daily_Data)-1)]
Beta_M_Daily = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Beta_M_Daily[,2:(dim(Daily_Data)[2]-1)] = as.matrix(Daily_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 252,min_obs = 90)[1,2,]/roll_cov(cbind(vwret,vwret),width = 252,min_obs = 90)[1,1,]),.SDcols=cols])

#Create matrix where we remove all down days for each company (set them to NA)
Daily_Data_Up = as.data.frame(Daily_Data)
setDT(Daily_Data_Up)
v1 <- colnames(Daily_Data_Up)[2:(ncol(Daily_Data_Up)-1)]
for(j in v1){
  set(Daily_Data_Up, i = which(Daily_Data_Up$vwret < 0), j = j, value = NA_integer_)
}
Beta_Up_Daily = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Beta_Up_Daily[,2:(dim(Daily_Data_Up)[2]-1)] = as.matrix(Daily_Data_Up[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 252,min_obs = 90)[1,2,]/roll_cov(cbind(vwret,vwret),width = 252,min_obs = 90)[1,1,]),.SDcols=cols])

#Create matrix where we remove all up days for each company (set them to NA)
Daily_Data_Down = as.data.frame(Daily_Data)
setDT(Daily_Data_Down)
v1 <- colnames(Daily_Data_Down)[2:(ncol(Daily_Data_Down)-1)]
for(j in v1){
  set(Daily_Data_Down, i = which(Daily_Data_Down$vwret > 0), j = j, value = NA_integer_)
}
Beta_Down_Daily = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Beta_Down_Daily[,2:(dim(Daily_Data_Down)[2]-1)] = as.matrix(Daily_Data_Down[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 252,min_obs = 90)[1,2,]/roll_cov(cbind(vwret,vwret),width = 252,min_obs = 90)[1,1,]),.SDcols=cols])



## Weighting:
h = log(2)/(2/3*252)
Weights = exp(-abs(252:1)*h)/(sum(exp(-abs(252:1)*h)))

Weighting_Beta_M_Daily = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_M_Daily[,2:(dim(Daily_Data)[2]-1)] = as.matrix(Daily_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 252,min_obs = 90, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 252,min_obs = 90, weights = Weights)[1,1,]),.SDcols=cols])

#Create matrix where we remove all down days for each company (set them to NA)
Daily_Data_Up = as.data.frame(Daily_Data)
setDT(Daily_Data_Up)
v1 <- colnames(Daily_Data_Up)[2:(ncol(Daily_Data_Up)-1)]
for(j in v1){
  set(Daily_Data_Up, i = which(Daily_Data_Up$vwret < 0), j = j, value = NA_integer_)
}
Weighting_Beta_Up_Daily = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_Up_Daily[,2:(dim(Daily_Data_Up)[2]-1)] = as.matrix(Daily_Data_Up[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 252,min_obs = 90, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 252,min_obs = 90, weights = Weights)[1,1,]),.SDcols=cols])

#Create matrix where we remove all down days for each company (set them to NA)
Daily_Data_Down = as.data.frame(Daily_Data)
setDT(Daily_Data_Down)
v1 <- colnames(Daily_Data_Down)[2:(ncol(Daily_Data_Down)-1)]
for(j in v1){
  set(Daily_Data_Down, i = which(Daily_Data_Down$vwret > 0), j = j, value = NA_integer_)
}
Weighting_Beta_Down_Daily = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_Down_Daily[,2:(dim(Daily_Data_Down)[2]-1)] = as.matrix(Daily_Data_Down[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 252,min_obs = 90, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 252,min_obs = 90, weights = Weights)[1,1,]),.SDcols=cols])


## Weighting - 60 Months:
h = log(2)/(2/3*1260)
Weights = exp(-abs(1260:1)*h)/(sum(exp(-abs(1260:1)*h)))

Weighting_Beta_M_Daily_60M = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_M_Daily_60M[,2:(dim(Daily_Data)[2]-1)] = as.matrix(Daily_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 1260,min_obs = 90, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 1260,min_obs = 90, weights = Weights)[1,1,]),.SDcols=cols])

#Create matrix where we remove all down days for each company (set them to NA)
Daily_Data_Up = as.data.frame(Daily_Data)
setDT(Daily_Data_Up)
v1 <- colnames(Daily_Data_Up)[2:(ncol(Daily_Data_Up)-1)]
for(j in v1){
  set(Daily_Data_Up, i = which(Daily_Data_Up$vwret < 0), j = j, value = NA_integer_)
}
Weighting_Beta_Up_Daily_60M = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_Up_Daily_60M[,2:(dim(Daily_Data_Up)[2]-1)] = as.matrix(Daily_Data_Up[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 1260,min_obs = 90, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 1260,min_obs = 90, weights = Weights)[1,1,]),.SDcols=cols])

#Create matrix where we remove all down days for each company (set them to NA)
Daily_Data_Down = as.data.frame(Daily_Data)
setDT(Daily_Data_Down)
v1 <- colnames(Daily_Data_Down)[2:(ncol(Daily_Data_Down)-1)]
for(j in v1){
  set(Daily_Data_Down, i = which(Daily_Data_Down$vwret > 0), j = j, value = NA_integer_)
}
Weighting_Beta_Down_Daily_60M = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_Down_Daily_60M[,2:(dim(Daily_Data_Down)[2]-1)] = as.matrix(Daily_Data_Down[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 1260,min_obs = 90, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 1260,min_obs = 90, weights = Weights)[1,1,]),.SDcols=cols])



## Weighting - 24 Months:
h = log(2)/(2/3*504)
Weights = exp(-abs(504:1)*h)/(sum(exp(-abs(504:1)*h)))

Weighting_Beta_M_Daily_24M = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_M_Daily_24M[,2:(dim(Daily_Data)[2]-1)] = as.matrix(Daily_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 504,min_obs = 90, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 504,min_obs = 90, weights = Weights)[1,1,]),.SDcols=cols])

#Create matrix where we remove all down days for each company (set them to NA)
Daily_Data_Up = as.data.frame(Daily_Data)
setDT(Daily_Data_Up)
v1 <- colnames(Daily_Data_Up)[2:(ncol(Daily_Data_Up)-1)]
for(j in v1){
  set(Daily_Data_Up, i = which(Daily_Data_Up$vwret < 0), j = j, value = NA_integer_)
}
Weighting_Beta_Up_Daily_24M = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_Up_Daily_24M[,2:(dim(Daily_Data_Up)[2]-1)] = as.matrix(Daily_Data_Up[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 504,min_obs = 90, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 504,min_obs = 90, weights = Weights)[1,1,]),.SDcols=cols])

#Create matrix where we remove all down days for each company (set them to NA)
Daily_Data_Down = as.data.frame(Daily_Data)
setDT(Daily_Data_Down)
v1 <- colnames(Daily_Data_Down)[2:(ncol(Daily_Data_Down)-1)]
for(j in v1){
  set(Daily_Data_Down, i = which(Daily_Data_Down$vwret > 0), j = j, value = NA_integer_)
}
Weighting_Beta_Down_Daily_24M = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_Down_Daily_24M[,2:(dim(Daily_Data_Down)[2]-1)] = as.matrix(Daily_Data_Down[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 504,min_obs = 90, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 504,min_obs = 90, weights = Weights)[1,1,]),.SDcols=cols])


## Weighting - 6 Months:
h = log(2)/(2/3*126)
Weights = exp(-abs(126:1)*h)/(sum(exp(-abs(126:1)*h)))

Weighting_Beta_M_Daily_6M = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_M_Daily_6M[,2:(dim(Daily_Data)[2]-1)] = as.matrix(Daily_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 126,min_obs = 45, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 126,min_obs = 45, weights = Weights)[1,1,]),.SDcols=cols])

#Create matrix where we remove all down days for each company (set them to NA)
Daily_Data_Up = as.data.frame(Daily_Data)
setDT(Daily_Data_Up)
v1 <- colnames(Daily_Data_Up)[2:(ncol(Daily_Data_Up)-1)]
for(j in v1){
  set(Daily_Data_Up, i = which(Daily_Data_Up$vwret < 0), j = j, value = NA_integer_)
}
Weighting_Beta_Up_Daily_6M = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_Up_Daily_6M[,2:(dim(Daily_Data_Up)[2]-1)] = as.matrix(Daily_Data_Up[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 126,min_obs = 45, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 126,min_obs = 45, weights = Weights)[1,1,]),.SDcols=cols])

#Create matrix where we remove all down days for each company (set them to NA)
Daily_Data_Down = as.data.frame(Daily_Data)
setDT(Daily_Data_Down)
v1 <- colnames(Daily_Data_Down)[2:(ncol(Daily_Data_Down)-1)]
for(j in v1){
  set(Daily_Data_Down, i = which(Daily_Data_Down$vwret > 0), j = j, value = NA_integer_)
}
Weighting_Beta_Down_Daily_6M = matrix(NA,nrow=dim(Daily_Data)[1],ncol=dim(Daily_Data)[2])
Weighting_Beta_Down_Daily_6M[,2:(dim(Daily_Data_Down)[2]-1)] = as.matrix(Daily_Data_Down[ , lapply(.SD, function(X) roll_cov(cbind(X,vwret),width = 126,min_obs = 45, weights = Weights)[1,2,]/roll_cov(cbind(vwret,vwret),width = 126,min_obs = 45, weights = Weights)[1,1,]),.SDcols=cols])





#############  so reduce the estimated dataset to monthly observations: ############

library(data.table)
library(lubridate)


Beta_M_Daily[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Beta_M_Daily = Beta_M_Daily[,-ncol(Beta_M_Daily)]
Beta_M_Daily = data.frame(Beta_M_Daily)
library(stringr)
colnames(Beta_M_Daily)[1:ncol(Beta_M_Daily)] = c("Date",str_remove_all(cols,pattern="X"))
Beta_M_Daily$Date = as.Date(as.character(Beta_M_Daily$Date),format="%Y%m%d")
#Pick the last date that month
Beta_M_Daily = setDT(Beta_M_Daily)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Beta_Up_Daily[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Beta_Up_Daily = Beta_Up_Daily[,-ncol(Beta_Up_Daily)]
Beta_Up_Daily = data.frame(Beta_Up_Daily)
colnames(Beta_Up_Daily)[1:ncol(Beta_Up_Daily)] = c("Date",str_remove_all(cols,pattern="X"))
Beta_Up_Daily$Date = as.Date(as.character(Beta_Up_Daily$Date),format="%Y%m%d")
#Pick the last date that month
Beta_Up_Daily = setDT(Beta_Up_Daily)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Beta_Down_Daily[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Beta_Down_Daily = Beta_Down_Daily[,-ncol(Beta_Down_Daily)]
Beta_Down_Daily = data.frame(Beta_Down_Daily)
colnames(Beta_Down_Daily)[1:ncol(Beta_Down_Daily)] = c("Date",str_remove_all(cols,pattern="X"))
Beta_Down_Daily$Date = as.Date(as.character(Beta_Down_Daily$Date),format="%Y%m%d")
#Pick the last date that month
Beta_Down_Daily = setDT(Beta_Down_Daily)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Weighting_Beta_M_Daily[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_M_Daily = Weighting_Beta_M_Daily[,-ncol(Weighting_Beta_M_Daily)]
Weighting_Beta_M_Daily = data.frame(Weighting_Beta_M_Daily)
colnames(Weighting_Beta_M_Daily)[1:ncol(Weighting_Beta_M_Daily)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_M_Daily$Date = as.Date(as.character(Weighting_Beta_M_Daily$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_M_Daily = setDT(Weighting_Beta_M_Daily)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Weighting_Beta_Up_Daily[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_Up_Daily = Weighting_Beta_Up_Daily[,-ncol(Weighting_Beta_Up_Daily)]
Weighting_Beta_Up_Daily = data.frame(Weighting_Beta_Up_Daily)
colnames(Weighting_Beta_Up_Daily)[1:ncol(Weighting_Beta_Up_Daily)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_Up_Daily$Date = as.Date(as.character(Weighting_Beta_Up_Daily$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_Up_Daily = setDT(Weighting_Beta_Up_Daily)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Weighting_Beta_Down_Daily[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_Down_Daily = Weighting_Beta_Down_Daily[,-ncol(Weighting_Beta_Down_Daily)]
Weighting_Beta_Down_Daily = data.frame(Weighting_Beta_Down_Daily)
colnames(Weighting_Beta_Down_Daily)[1:ncol(Weighting_Beta_Down_Daily)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_Down_Daily$Date = as.Date(as.character(Weighting_Beta_Down_Daily$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_Down_Daily = setDT(Weighting_Beta_Down_Daily)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

#Different Weighting Methods:

Weighting_Beta_M_Daily_60M[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_M_Daily_60M = Weighting_Beta_M_Daily_60M[,-ncol(Weighting_Beta_M_Daily_60M)]
Weighting_Beta_M_Daily_60M = data.frame(Weighting_Beta_M_Daily_60M)
colnames(Weighting_Beta_M_Daily_60M)[1:ncol(Weighting_Beta_M_Daily_60M)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_M_Daily_60M$Date = as.Date(as.character(Weighting_Beta_M_Daily_60M$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_M_Daily_60M = setDT(Weighting_Beta_M_Daily_60M)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Weighting_Beta_Up_Daily_60M[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_Up_Daily_60M = Weighting_Beta_Up_Daily_60M[,-ncol(Weighting_Beta_Up_Daily_60M)]
Weighting_Beta_Up_Daily_60M = data.frame(Weighting_Beta_Up_Daily_60M)
colnames(Weighting_Beta_Up_Daily_60M)[1:ncol(Weighting_Beta_Up_Daily_60M)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_Up_Daily_60M$Date = as.Date(as.character(Weighting_Beta_Up_Daily_60M$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_Up_Daily_60M = setDT(Weighting_Beta_Up_Daily_60M)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Weighting_Beta_Down_Daily_60M[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_Down_Daily_60M = Weighting_Beta_Down_Daily_60M[,-ncol(Weighting_Beta_Down_Daily_60M)]
Weighting_Beta_Down_Daily_60M = data.frame(Weighting_Beta_Down_Daily_60M)
colnames(Weighting_Beta_Down_Daily_60M)[1:ncol(Weighting_Beta_Down_Daily_60M)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_Down_Daily_60M$Date = as.Date(as.character(Weighting_Beta_Down_Daily_60M$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_Down_Daily_60M = setDT(Weighting_Beta_Down_Daily_60M)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 


Weighting_Beta_M_Daily_24M[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_M_Daily_24M = Weighting_Beta_M_Daily_24M[,-ncol(Weighting_Beta_M_Daily_24M)]
Weighting_Beta_M_Daily_24M = data.frame(Weighting_Beta_M_Daily_24M)
colnames(Weighting_Beta_M_Daily_24M)[1:ncol(Weighting_Beta_M_Daily_24M)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_M_Daily_24M$Date = as.Date(as.character(Weighting_Beta_M_Daily_24M$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_M_Daily_24M = setDT(Weighting_Beta_M_Daily_24M)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Weighting_Beta_Up_Daily_24M[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_Up_Daily_24M = Weighting_Beta_Up_Daily_24M[,-ncol(Weighting_Beta_Up_Daily_24M)]
Weighting_Beta_Up_Daily_24M = data.frame(Weighting_Beta_Up_Daily_24M)
colnames(Weighting_Beta_Up_Daily_24M)[1:ncol(Weighting_Beta_Up_Daily_24M)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_Up_Daily_24M$Date = as.Date(as.character(Weighting_Beta_Up_Daily_24M$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_Up_Daily_24M = setDT(Weighting_Beta_Up_Daily_24M)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Weighting_Beta_Down_Daily_24M[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_Down_Daily_24M = Weighting_Beta_Down_Daily_24M[,-ncol(Weighting_Beta_Down_Daily_24M)]
Weighting_Beta_Down_Daily_24M = data.frame(Weighting_Beta_Down_Daily_24M)
colnames(Weighting_Beta_Down_Daily_24M)[1:ncol(Weighting_Beta_Down_Daily_24M)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_Down_Daily_24M$Date = as.Date(as.character(Weighting_Beta_Down_Daily_24M$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_Down_Daily_24M = setDT(Weighting_Beta_Down_Daily_24M)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 



Weighting_Beta_M_Daily_6M[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_M_Daily_6M = Weighting_Beta_M_Daily_6M[,-ncol(Weighting_Beta_M_Daily_6M)]
Weighting_Beta_M_Daily_6M = data.frame(Weighting_Beta_M_Daily_6M)
colnames(Weighting_Beta_M_Daily_6M)[1:ncol(Weighting_Beta_M_Daily_6M)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_M_Daily_6M$Date = as.Date(as.character(Weighting_Beta_M_Daily_6M$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_M_Daily_6M = setDT(Weighting_Beta_M_Daily_6M)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Weighting_Beta_Up_Daily_6M[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_Up_Daily_6M = Weighting_Beta_Up_Daily_6M[,-ncol(Weighting_Beta_Up_Daily_6M)]
Weighting_Beta_Up_Daily_6M = data.frame(Weighting_Beta_Up_Daily_6M)
colnames(Weighting_Beta_Up_Daily_6M)[1:ncol(Weighting_Beta_Up_Daily_6M)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_Up_Daily_6M$Date = as.Date(as.character(Weighting_Beta_Up_Daily_6M$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_Up_Daily_6M = setDT(Weighting_Beta_Up_Daily_6M)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 

Weighting_Beta_Down_Daily_6M[,1] = Daily_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_Down_Daily_6M = Weighting_Beta_Down_Daily_6M[,-ncol(Weighting_Beta_Down_Daily_6M)]
Weighting_Beta_Down_Daily_6M = data.frame(Weighting_Beta_Down_Daily_6M)
colnames(Weighting_Beta_Down_Daily_6M)[1:ncol(Weighting_Beta_Down_Daily_6M)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_Down_Daily_6M$Date = as.Date(as.character(Weighting_Beta_Down_Daily_6M$Date),format="%Y%m%d")
#Pick the last date that month
Weighting_Beta_Down_Daily_6M = setDT(Weighting_Beta_Down_Daily_6M)[order(Date), .SD[which.max(Date)], by = .(format(Date,"%Y%m"))] 


save(Weighting_Beta_M_Daily_60M,file = "Weighting_Daily_Beta_M_60M.RData")
save(Weighting_Beta_Up_Daily_60M,file = "Weighting_Daily_Beta_Up_60M.RData")
save(Weighting_Beta_Down_Daily_60M,file = "Weighting_Daily_Beta_Down_60M.RData")

save(Weighting_Beta_M_Daily_24M,file = "Weighting_Daily_Beta_M_24M.RData")
save(Weighting_Beta_Up_Daily_24M,file = "Weighting_Daily_Beta_Up_24M.RData")
save(Weighting_Beta_Down_Daily_24M,file = "Weighting_Daily_Beta_Down_24M.RData")

save(Weighting_Beta_M_Daily_6M,file = "Weighting_Daily_Beta_M_6M.RData")
save(Weighting_Beta_Up_Daily_6M,file = "Weighting_Daily_Beta_Up_6M.RData")
save(Weighting_Beta_Down_Daily_6M,file = "Weighting_Daily_Beta_Down_6M.RData")


###################################################### Shrinkage #############################################################################


#Load in the prior results:

#Use as prior the cross-sectional average beta at that point in time:
load("Daily_Beta_M.RData")
load("Daily_Beta_Up.RData")
load("Daily_Beta_Down.RData")
#This is the Cross-section
Shrinking_Beta_M_Daily = matrix(NA,nrow=dim(Beta_M_Daily)[1],ncol=dim(Beta_M_Daily)[2])
Shrinking_Beta_Up_Daily = matrix(NA,nrow=dim(Beta_M_Daily)[1],ncol=dim(Beta_M_Daily)[2])
Shrinking_Beta_Down_Daily = matrix(NA,nrow=dim(Beta_M_Daily)[1],ncol=dim(Beta_M_Daily)[2])
library(data.table)
setDF(Beta_M_Daily)
setDF(Beta_Up_Daily)
setDF(Beta_Down_Daily)
#For Beta_M
cross_Beta_M = rowMeans(Beta_M_Daily[,2:ncol(Beta_M_Daily)], na.rm=T)
sd_cross_Beta_M = (sd(rowMeans(Beta_M_Daily[,2:ncol(Beta_M_Daily)], na.rm=T),na.rm=T))^2
prior_Beta_M = (apply(Beta_M_Daily[,2:ncol(Beta_M_Daily)],1,sd,na.rm=T))^2
#For Beta_Up
cross_Beta_Up = rowMeans(Beta_Up_Daily[,2:ncol(Beta_Up_Daily)], na.rm=T)
sd_cross_Beta_Up = (sd(rowMeans(Beta_Up_Daily[,2:ncol(Beta_Up_Daily)], na.rm=T),na.rm=T))^2
prior_Beta_Up = (apply(Beta_Up_Daily[,2:ncol(Beta_Up_Daily)],1,sd,na.rm=T))^2
#For Beta_Down
cross_Beta_Down = rowMeans(Beta_Down_Daily[,2:ncol(Beta_Down_Daily)], na.rm=T)
sd_cross_Beta_Down = (sd(rowMeans(Beta_Down_Daily[,2:ncol(Beta_Down_Daily)], na.rm=T),na.rm=T))^2
prior_Beta_Down = (apply(Beta_Down_Daily[,2:ncol(Beta_Down_Daily)],1,sd,na.rm=T))^2

#Do we shrink this based on the entire sample, apparently (...)
for(i in 2:ncol(Beta_M_Daily)){
  Shrinking_Beta_M_Daily[,i] = (Beta_M_Daily[,i] + sd_cross_Beta_M/prior_Beta_M * cross_Beta_M)/(1 + sd_cross_Beta_M/prior_Beta_M)

  Shrinking_Beta_Up_Daily[,i] = (Beta_Up_Daily[,i] + sd_cross_Beta_Up/prior_Beta_Up * cross_Beta_Up)/(1 + sd_cross_Beta_Up/prior_Beta_Up)
  
  Shrinking_Beta_Down_Daily[,i] = (Beta_Down_Daily[,i] + sd_cross_Beta_Down/prior_Beta_Down * cross_Beta_Down)/(1 + sd_cross_Beta_Down/prior_Beta_Down)
  
  print(paste0(i," of iteration ",ncol(Beta_M_Daily)," which is ",i/ncol(Beta_M_Daily)*100," Percent"))
}
#Add the Dates:
Shrinking_Beta_M_Daily = data.frame(Shrinking_Beta_M_Daily)
Shrinking_Beta_M_Daily[,1] = Beta_M_Daily[,1]
colnames(Shrinking_Beta_M_Daily) = colnames(Beta_M_Daily)

Shrinking_Beta_Up_Daily = data.frame(Shrinking_Beta_Up_Daily)
Shrinking_Beta_Up_Daily[,1] = Beta_Up_Daily[,1]
colnames(Shrinking_Beta_Up_Daily) = colnames(Beta_Up_Daily)

Shrinking_Beta_Down_Daily = data.frame(Shrinking_Beta_Down_Daily)
Shrinking_Beta_Down_Daily[,1] = Beta_Down_Daily[,1]
colnames(Shrinking_Beta_Down_Daily) = colnames(Beta_M_Daily)



load("Weighting_Daily_Beta_M.RData")
load("Weighting_Daily_Beta_Up.RData")
load("Weighting_Daily_Beta_Down.RData")
#This is the Cross-section
Weighting_Shrinking_Beta_M_Daily = matrix(NA,nrow=dim(Weighting_Beta_M_Daily)[1],ncol=dim(Weighting_Beta_M_Daily)[2])
Weighting_Shrinking_Beta_Up_Daily = matrix(NA,nrow=dim(Weighting_Beta_M_Daily)[1],ncol=dim(Weighting_Beta_M_Daily)[2])
Weighting_Shrinking_Beta_Down_Daily = matrix(NA,nrow=dim(Weighting_Beta_M_Daily)[1],ncol=dim(Weighting_Beta_M_Daily)[2])
library(data.table)
setDF(Weighting_Beta_M_Daily)
setDF(Weighting_Beta_Up_Daily)
setDF(Weighting_Beta_Down_Daily)

#For Beta_M
cross_Beta_M = rowMeans(Weighting_Beta_M_Daily[,2:ncol(Weighting_Beta_M_Daily)], na.rm=T)
sd_cross_Beta_M = (sd(rowMeans(Weighting_Beta_M_Daily[,2:ncol(Weighting_Beta_M_Daily)], na.rm=T),na.rm=T))^2
prior_Beta_M = (apply(Weighting_Beta_M_Daily[,2:ncol(Weighting_Beta_M_Daily)],1,sd,na.rm=T))^2
#For Beta_Up
cross_Beta_Up = rowMeans(Weighting_Beta_Up_Daily[,2:ncol(Weighting_Beta_Up_Daily)], na.rm=T)
sd_cross_Beta_Up = (sd(rowMeans(Weighting_Beta_Up_Daily[,2:ncol(Weighting_Beta_Up_Daily)], na.rm=T),na.rm=T))^2
prior_Beta_Up = (apply(Weighting_Beta_Up_Daily[,2:ncol(Weighting_Beta_Up_Daily)],1,sd,na.rm=T))^2
#For Beta_Down
cross_Beta_Down = rowMeans(Weighting_Beta_Down_Daily[,2:ncol(Weighting_Beta_Down_Daily)], na.rm=T)
sd_cross_Beta_Down = (sd(rowMeans(Weighting_Beta_Down_Daily[,2:ncol(Weighting_Beta_Down_Daily)], na.rm=T),na.rm=T))^2
prior_Beta_Down = (apply(Weighting_Beta_Down_Daily[,2:ncol(Weighting_Beta_Down_Daily)],1,sd,na.rm=T))^2

for(i in 2:ncol(Beta_M_Daily)){
  Weighting_Shrinking_Beta_M_Daily[,i] = (Weighting_Beta_M_Daily[,i] + sd_cross_Beta_M/prior_Beta_M * cross_Beta_M)/(1 + sd_cross_Beta_M/prior_Beta_M)
  
  Weighting_Shrinking_Beta_Up_Daily[,i] = (Weighting_Beta_Up_Daily[,i] + sd_cross_Beta_Up/prior_Beta_Up * cross_Beta_Up)/(1 + sd_cross_Beta_Up/prior_Beta_Up)
  
  Weighting_Shrinking_Beta_Down_Daily[,i] = (Weighting_Beta_Down_Daily[,i] + sd_cross_Beta_Down/prior_Beta_Down * cross_Beta_Down)/(1 + sd_cross_Beta_Down/prior_Beta_Down)
  
  print(paste0(i," of iteration ",ncol(Weighting_Beta_M_Daily)," which is ",i/ncol(Weighting_Beta_M_Daily)*100," Percent"))
}
#Add the Dates:
Weighting_Shrinking_Beta_M_Daily = data.frame(Weighting_Shrinking_Beta_M_Daily)
Weighting_Shrinking_Beta_M_Daily[,1] = Weighting_Beta_M_Daily[,1]
colnames(Weighting_Shrinking_Beta_M_Daily) = colnames(Weighting_Beta_M_Daily)

Weighting_Shrinking_Beta_Up_Daily = data.frame(Weighting_Shrinking_Beta_Up_Daily)
Weighting_Shrinking_Beta_Up_Daily[,1] = Weighting_Beta_Up_Daily[,1]
colnames(Weighting_Shrinking_Beta_Up_Daily) = colnames(Weighting_Beta_Up_Daily)

Weighting_Shrinking_Beta_Down_Daily = data.frame(Weighting_Shrinking_Beta_Down_Daily)
Weighting_Shrinking_Beta_Down_Daily[,1] = Weighting_Beta_Down_Daily[,1]
colnames(Weighting_Shrinking_Beta_Down_Daily) = colnames(Weighting_Beta_Down_Daily)



############################################### Save the Data #########################################
save(Beta_M_Daily,file = "Daily_Beta_M.RData")
save(Beta_Up_Daily,file = "Daily_Beta_Up.RData")
save(Beta_Down_Daily,file = "Daily_Beta_Down.RData")

save(Weighting_Beta_M_Daily,file = "Weighting_Daily_Beta_M.RData")
save(Weighting_Beta_Up_Daily,file = "Weighting_Daily_Beta_Up.RData")
save(Weighting_Beta_Down_Daily,file = "Weighting_Daily_Beta_Down.RData")

save(Shrinking_Beta_M_Daily, file = "Shrinking_Beta_M_Daily.RData")
save(Shrinking_Beta_Up_Daily, file = "Shrinking_Beta_Up_Daily.RData")
save(Shrinking_Beta_Down_Daily, file = "Shrinking_Beta_Down_Daily.RData")


save(Weighting_Shrinking_Beta_M_Daily, file = "Weighting_Shrinking_Beta_M_Daily.RData")
save(Weighting_Shrinking_Beta_Up_Daily, file = "Weighting_Shrinking_Beta_Up_Daily.RData")
save(Weighting_Shrinking_Beta_Down_Daily, file = "Weighting_Shrinking_Beta_Down_Daily.RData")



